Slow dynamics and precursors of the glass transition in granular fluids 
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We use event driven simulations to analyze glassy dynamics as a function of density and energy 
dissipation in a two-dimensional bidisperse granular fluid under stationary conditions. Clear signa- 
tures of a glass transition are identified, such as an increase of relaxation times over several orders 
of magnitude. As the inelasticity is increased, the glass transition is shifted to higher densities 
and the precursors of the transition become less and less pronounced - in agreement with a recent 
mode-coupling theory. We analyze the long-time tails of the velocity autocorrelation and discuss its 
consequences for the nonexistence of the diffusion constant in two dimensions. 
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I. INTRODUCTION 

The physics of nonequilibrium systems in general in- 
volves dissipation and energy injection. These proper- 
ties distinguish nonequilibrium but stationary systems 
from those in thermodynamic equilibrium and render the 
application of standard thermodynamics and statistical 
mechanics impossible such that one has to develop new 
methods to treat nonequilibrium systems or resort to ex- 
perimental and computational methods. Due to their 
inherent dissipative nature, driven granular fluids serve 
as model systems for nonequilibrium phenomena and in 
particular phase transitions in macroscopic systems out 
of equilibrium. Questions addressed include the crystal- 
lization of granular fluids [1, 2], the jamming transition 
[3] and the glass transition [4] . 

Here we focus on a two-dimensional driven granular 
fluid, displaying structural arrest as the density is in- 
creased. Most experimental work in the field does indeed 
use a two-dimensional setup, allowing for extensive parti- 
cle tracking with high speed and resolution cameras. Sev- 
eral driving mechanisms have been studied, such as ver- 
tical vibrations [1, 2, 5] and air-tables [6-8]. Since we are 
interested in structural arrest, the most relevant experi- 
ments are those of Abate and Durian [6, 7], who showed 
that the approach to jamming in a granular medium re- 
sembles the glass transition in a molecular or colloidal 
fluid: The mean square displacement develops a plateau 
and the dynamics becomes strongly heterogeneous. 

Recently a mode-coupling theory [4] has been set up 
for a driven granular fluid in the stationary state. A 
glass transition was located below random close packing 
which is qualitatively similar, but quantitatively different 
from the corresponding transition for a molecular fluid. 
The coherent scattering function was shown to display a 
plateau - provided the density is sufficiently close to the 
glass transition. The timescale for relaxation from the 



plateau (a-relaxation) diverges at the glass transition. 
The transition density as well as the dynamics depend on 
the degree of inelasticity, parametrized by the coefficient 
of restitution. 

In this paper we investigate the dynamics of a driven 
two-dimensional granular fluid close to the glass transi- 
tion, mimicking an air-table experiment. We use event 
driven simulations to compute the mean square displace- 
ment, the velocity autocorrelation, and the incoherent 
scattering function. 



II. MODEL 

Dense systems of monodisperse, elastic hard disks are 
known to undergo a phase transition to a crystalline state 
[9]. The only state variable is the packing fraction r\, 
defined by the area covered by particles divided by the 
area of the whole system. The phase transition occurs for 
packing fractions between 77/ = 0.691 and r\ s = 0.7163, 
where packing fractions up to 77/ are completely in the 
fluid phase and packing fractions larger than r) s are com- 
pletely in the solid phase. Packing fractions between 
r/f and rj s fall into the coexistence area. An overview 
over simulations performed to determine the phase tran- 
sition is given in [10]. Crystallization can be avoided 
in binary mixtures of particles with sufficiently differ- 
ent sizes. Whereas bidisperse mixtures of elastic hard 
spheres form stable crystalline structures for a size ratio 
smaller than 1.2 [11], larger size ratios suffice to avoid 
crystallization. We expect similar behaviour for a gran- 
ular fluid and choose a bidisperse symmetric mixture of 
hard disks. The ratio of the radius of the big particles, 
Rb, to the radius of the small particles, R s , is given by 
Rb/Rs = 1.43. 



Collision rules 
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We restrict our analysis to smooth particles without 
be addressed. rotational degree of freedom. Particles' positions and ve- 
locities are denoted by {n, Vi}£L x . The dissipation arises 
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solely from incomplete normal restitution with a constant 
coefficient of restitution e: The change of the relative ve- 
locity g := vi — V2 in the direction of the normal vector 
n := (n - r 2 )/|ri - r 2 | is given by 



fo- 



ri 



-e(g-n). 



(1) 



Here primed and unprimed quantities indicate postcol- 
lisional and precollisional velocities of the two colliding 
particles. The postcollisional velocities of the two col- 
liding disks are expressed in terms of the precollisional 
velocities according to 
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where 
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Here m, denotes the mass of particle i and we assume 
constant mass density for all particles such that the 
mass ratio of the big to the small particles is given by 
(Rb/R s ) 2 ■ 



B. Driving 

Due to inelastic collisions, energy is dissipated and we 
have to feed energy into the system in order to obtain a 
stationary state. One of the simplest driving mechanisms 
[12] is to kick a given particle, say, particle i, instanta- 
neously at time t according to 



(5) 



Here pr> r is the driving amplitude and (t) indicates the 
direction of the driving which is chosen randomly with 
{€i(t)tj(t')) = 5 t] 5 af} 5{t - t'). The Cartesian compo- 
nents £f,a = x,y, are distributed according to a Gaus- 
sian distribution with zero mean. 

In practice, we implement the stochastic driving pro- 
cess by kicking the particles randomly with frequency 
for- In order to conserve momentum locally, we apply 
kicks of equal magnitude and opposite direction to pairs 
of neighbouring particles [13]. Thereby momentum is 
conserved locally. 



We use an efficient standard event driven algorithm 
[14] that enables us to simulate large systems for long 
times. Inelastically colliding hard particles are known to 
undergo the so-called inelastic collapse. We circumvent 
the inelastic collapse with the same method used in [15]. 

We consider a system of N = 10000 particles in an 
area A = L 2 , where L denotes the length of system. 
The mixture is chosen symmetric, so that small and big 
particles appear in equal numbers. Initial velocities are 
drawn from a Gaussian distribution. To analyze the 
slow dynamics of the system as the glass transition is 
approached, we simulate a range of coefficients of resti- 
tution, 0.5 < e < 0.95 and area fractions 0.1 < r\ < 0.81. 
For every parameter set (77, e) we performed an event 
driven molecular dynamic simulation with approximately 
10 7 collisions per particle. 

Several initial configurations with r] = 0.84 were taken 
from a time driven molecular dynamics simulation of soft 
disks [16], that suffer dissipation and finally come to rest. 
By expanding the system we change the packing fraction, 
allowing us to investigate systems of hard disks up to 
rj = 0.81. To be sure that our simulations do not depend 
on the initial configurations (i.e., the initial values of 
particles' velocities and positions), we have repeated the 
simulation for each parameter set, using several different 
initial configurations. 

In order to investigate the long-time behaviour of the 
mean square displacement and the velocity autocorrela- 
tion function it was necessary to simulate larger systems 
with N = 4 • 10 6 . This was done for e = 0.7, 0.8, 0.9 and 
r] = 0.1... 0.78. 



III. RESULTS 

Our main concern is the slow dynamics in a dense gran- 
ular fluid and in particular the precursors of the glass 
transition. Our analysis will be based on the time depen- 
dent van Hove correlation function, the mean square dis- 
placemnet (MSD) and the velocity autocorrelation func- 
tion (VACF) in the stationary state. Hence we briefly 
regress to discuss the relaxation towards stationarity. 



C. Implementation 

As long as the system is not jammed, the collisions 
of hard particles are instantaneous, so that event driven 
simulations can be applied. Lasting contacts, that we ex- 
pect to occur in a jammed state, cannot be accounted for 
by event driven simulations, so that we are restricted to 
densities below the jamming point. The glass transition 
in a driven granular fluid is expected to occur at den- 
sities below jamming, so that the restriction is actually 
less severe. 



A. Stationary state 

To maintain a stationary state, the dissipated energy 
due to inelastic collisions must be compensated by a com- 
parable driving energy fed into the system. In a two- 
dimensional monodisperse system, we define the granular 
temperature as the average kinetic energy 

1 N 

T := — > — v 2 (6) 
N ^ 2 ' 

i=l 
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A stationary state can be achieved by balancing dissipa- 
tion and driving [17] 
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The collisional loss is estimated by the energy loss in a 
i — 2 

single collision, 4 £ T, multiplied by the Enskog collision 
frequency uj e : 



oj e = 2Trg(2a)(2a) d ~ 1 



NJTU) 



(8) 



Here g (2a) is the pair correlation at contact. 

In the following we are going to discuss correlations, 
such as the MSD or the VACF as a function of density 
and inelasticity e. We are particularly interested in the 
approach to the glass transition as the inelasticity is var- 
ied. To isolate effects of inelasticity and packing fraction, 
we have to keep the temperature constant (without loss of 
generality we take T ~ 1). This can be achieved approx- 
imately by adjusting the driving (see Eq. 7) according 
to 
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In order to sustain a stationary state we choose the driv- 
ing frequency of the order of the collision frequency as 
suggested by Bizon et al. [18]. Note that the driving am- 
plitude scales with the coefficient of restitution, such that 
we are able to take the elastic limit, e — >■ 1. 

In a symmetric binary mixture, we define two temper- 
atures, one for the big particles (Tj) and one for the small 
ones (T s ): 



n 

T, 



(10) 



These are in general unequal in the stationary state due 
to violation of equipartition [19]. The analogous cool- 
ing equations for and T s have been derived, also in 
the presence of driving, and can in principle be used to 
determine the driving amplitude [20]. We refrain from 
an exact solution and instead choose the driving ampli- 
tude por ~ Vl — £ 2 and absorb the density dependence 
in /j> r which is chosen equal to the collision frequency. 
An example for the relaxation to the stationary state is 
shown in Fig.l, where we plot the two temperatures Tb 
and T s . 

In the following we use units of length, mass and time, 
such that the radius of the small particles is one, their 
mass is one and |mjV? (i = 0) = 1 (average over all parti- 
cles). The natural time step in an event driven simulation 
is the number of collisions. In the stationary state the 




FIG. 1: (Color online) Time dependence of the temperature 
of the small particles (T s , upper curve) and big particles (Tj, 
lower curve) for e = 0.90 and packing fraction 77 = 0.5. Inset: 
Number of collisions per particle s as a function of time in 
dimensionless units. 



number of collisions increases linearly in time, as shown 
in the inset of Fig. 1 so that the two are related by a 
constant factor. 



B. Long-time tails of the Velocity Autocorrelation 



The VACF is defined by 



1 N 



(ii) 



where the average (...) is taken over initial conditions. 
The VACF is known to exhibit a hydrodynamic long-time 
tail in systems of hard spheres [21, 22] and hard disks [23] . 
Such long-time tails have also been observed in driven 
granular fluids in three dimensions. It was furthermore 
shown by simulations and a simple scaling argument that 
the tails depend sensitively on the driving mechanism: 
If driving conserves momentum locally, then the i~ 3 / 2 - 
tail known from elastic systems prevails, whereas it is 
replaced by t , if momentum is not conserved by the 
driving [15]. 

For two-dimensional fluids the situation is less clear. 
Two functional forms of the long-time decay of the 
equilibrium VACF have been proposed. Hydrodynamic 
theory [24] as well as mode-coupling theory [25] pre- 
dict $hd oc t -1 . If these approaches are made self- 
consistent [26, 27], the time decay is changed to <!>sc oc 
i~ 1 (y / m"(i)) _1 in agreement with renormalization group 
theory [28]. Recent simulations [23] claim that the self 
consistent theory is supported by large scale simulations. 

We compute the VACF for a bidisperse mixture of 
N = 4 • 10 6 particles as suggested by [23] in order to 
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FIG. 2: (Color online) The VACF for small particles for low 
and intermediate packing fractions as a function of time mea- 
sured in number of collisions per particle s. The long-time 
tails are more pronounced for increasing r\. 



FIG. 3: (Color online) The VACF (lower curves) for small par- 
ticles for r\ = 0.4 (left) and r\ = 0.5 (right) and e = 0.7, divided 
by both predictions $ oc s _1 (middle), oc s _ (yhi(i)) _ (up- 
per curves); time given in number of collisions per particle 
s. 



be able to analyze the long-time behavior of the VACF. 
The more inelastic system, e — 0.7, is expected to ex- 
hibit the most prominent long-time tail as observed in 
three-dimensional inelastic fluids [15]. Fig. 2 shows the 
VACF for an inelastic mixture with e = 0.7 and packing 
fractions up to r\ — 0.5. The long-time tail in the decay 
seems to be most pronounced for 77 = 0.4 and 77 = 0.5. In 
order to test the above predictions, we divide the VACF 
by each of the suggested decay laws and show the result 
in Fig. 3. If one of the two decay laws described the decay 
exactly, then the resulting curve should fluctuate around 
a constant. This is not the case neither for oc t" 1 nor for 
cx £ _1 -y/ln(t) . Hence, neither of the suggestions is sup- 
ported by our data, however we cannot exclude that one 
has to go to even longer times to identify the long-time 
tail unambiguously. 

For increasing e and high packing fractions the VACF 
is expected to show backscattering effects, i.e., a time 
interval, in which the VACF becomes negative because 
the particles are reflected from the cage formed by their 
neighbours. This is confirmed by the VACF in a system 
of 77 = 0.72 where the negative part of the VACF sets 
in for larger e after showing a dip for intermediate s, see 
Fig. 4. For constant e, the range of the negative VACF 
is increased for increasing 77, see inset of Fig. 4. 

C. Mean Square Displacement and Diffusion 
Coefficients 

One prominent signature of the glass transition is the 
formation of a plateau in the MSD 

1 N 

((Ar(t)) 2 ) = -^((r 4 (i)-r 4 (0)) 2 ). (12) 
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FIG. 4: (Color online) The VACF for small particles for dif- 
ferent e at packing fraction 77 = 0.72. Marked data points 
indicate a negative value. Inset: The VACF for e = 0.7 and 
packing fractions r\ = 0.76, 0.78; time given in number of col- 
lisions per particle s. 



The MSD for a mildly inelastic system (coefficient of 
restitution e = 0.9) is shown in Fig. 5 for a wide range of 
packing fractions 0.1 < 77 < 0.81. The simulated system 
consists of 10000 particles, which have been simulated for 
10 7 collisions per particle. Subsequently, the data was 
divided into approximately 1000 time intervals each of 
them giving rise to a single MSD. These were averaged to 
obtain the final time dependent MSD [29] with a standard 
error of 3 ■ 10~ 4 . Ballistic motion is observed for small 
times for all densities. The larger the density the smaller 
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FIG. 5: (Color online) MSD for small particles for e = 0.90 
and different packing fractions r\ — 0.1; 0.2; 0.3; 0.4; 0.5; 0.6; 
0.65; 0.7; 0.72; 0.74; 0.75; 0.76; 0.77; 0.78; 0.79; 0.80; 0.81. 
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FIG. 6: (Color online) MSD for small particles for e = 0.70 
and different packing fractions r\ — 0.1; 0.2; 0.3; 0.4; 0.5; 0.6; 
0.65; 0.7; 0.72; 0.74; 0.75; 0.76; 0.77; 0.78; 0.79; 0.80; 0.81. 



the ballistic regime. For packing fractions rj < 0.7, bal- 
listic motion crosses over to diffusive behaviour at long 
times. For rj > 0.7, a plateau develops, indicating that a 
cage has formed which restricts the particles motion to 
the space within the cage. Actually, our data indicate 
subdiffusive behavior instead of a plateau; a blowup of 
the relevant time range is shown in the inset of Fig. 5. 
For times 1 < t < 10 the MSDs for 77 = 0.81 show alge- 
braic increase with exponents 0.6163(11) for e = 0.7 and 
0.547(6) for e = 0.9, respectively. 

For area fraction 0.7 < 77 < 0.77 the cage is observed 
to break up for longer times so that ultimately diffusion 
prevails except for the highest density. 

As the system becomes more inelastic we expect the 
glass transition to shift to higher densities [4]. Further- 
more the range of densities where a clear separation of 
timescales can be observed, is expected to shrink. These 
expectations are indeed born out by our simulations. The 
MSD for a more inelastic system (coefficient of restitu- 
tion e = 0.7) is shown in Fig. 6. The plateau is hardly 
visible as compared to the more elastic system, instead 
subdiffusive behaviour is observed for the highest densi- 
ties. 

Next, we want to quantify the prediction of granular 
mode-coupling theory (GMCT) that the more inelastic 
systems require increasing critical densities for the glass 
transition to occur. This is done in two steps. We first 
identify a timescale, such that the small particles' MSD 
= 1. This corresponds to the time at which the mean 
travelled distance of all small particles equals their ra- 
dius. This timescale, i.e. t(MSD = 1), is shown in Fig. 
7 as a function of area fraction for a wide range of e. 
As can be seen in the figure, this time scale strongly in- 
creases as the glass transition is approached, the more 
so the more elastic the system. Now we ask: at what 
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FIG. 7: (Color online) Times t(r], e) of equal MSD for small 
particles for s = 0.5, 0.55, 0.6, 0.65, 0.7, 0.75, 0.8, 0.85, 0.9, 
0.95. Inset: Lines of equal MSDs for different times. 



area fraction has the timescale acquired a given value, 
e.g. t(MSD = 1) = 3? This yields an area fraction as 
a function of e, which is shown as the blue curve in the 
inset and corresponds to the lowest broken line in the 
main part of the figure. Hence, in the inset we show lines 
of equal relaxation time (t = 3,5,8) in the (77, e) plane. 
For all times t the packing fraction 77 is a monotonically 
decreasing function of e in agreement with GMCT [4] . 
The MSD is related to the VACF by ^Ar 2 (t) = 

2 Jj $(t')dt'. With this relat ion we can define a time de- 
pendent diffusion coefficient D(t) := 1/2 J* $(i')di'. The 
predicted functional forms of the VACF, $hd (t) oc t^ 1 
and 3>sc(i) oc t -1 (^/ln(f)) -1 , respectively, yield the fol- 
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FIG. 8: (Color online) Time dependent diffusion coefficient 
D(s) for small particles for e — 0.7 and different packing 
fractions. 
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FIG. 9: (Color online) Full lines: D(s)/ln(s), dashed lines: 



lowing time dependent diffusion coefficients: 

DhdW oc lnt, (13) 
Dsc(t) cx Vml. (14) 

Both of these expressions diverge for t — > oo, hence for 
both the long-time limit of the diffusion coefficient does 
not exist. 

We first check that D(s) = ±^(Ar(s)) 2 ) does not at- 
tain a stationary value but increases indefinitely as a 
function of time. Data for e = 0.7 are shown in Fig. 
8 for a wide range of packing fractions. Obviously no 
stationary value is attained for any packing fraction; the 
same holds for other values of e. 

In order to check if one of the proposed time depen- 
dences in Eqs. (13) or (14) is correct, we plot D(s)/ m(s) 
and D(s)/ ^ln(s). Similar to the test of the VACF, these 
quantities should fluctuate around a constant if they were 
describing D(s) correctly. The resulting graphs are plot- 
ted in Fig. 9. 

Neither the prediction of MCT, D(t) cx y/ln(t), nor 
the prediction of hydrodynamics, D{t) cx ln(i), account 
for the data in the full range of densities. MCT works 
slightly better at higher densities, whereas hydrodynam- 
ics fits the data better at intermediate densities. Hence, 
our data are not sufficient to discriminate between the 
two alternatives. 

The above discussion shows that the long-time behav- 
ior in a two-dimensional fluid is not truely diffusive. On 
the other hand it is known that the long-time tails of the 
VACF are suppressed close to the glass transition, where 
the shear viscosity v increases dramatically. In fact, the 
MCT of Ernst et al. [25] predicts for a two-dimensional 
clastic fluid 

kT 1 

®(t)oc- 7 — -. 15 

v ; 8rmr(i/ + D) t ' 



For large packing fractions in the limit of the glass transi- 
tion, v diverges and thus, the long-time tail in the VACF 
is suppressed. This allows us to calculate the diffusion 
coefficient approximately by two methods: 

• We calculate the time dependent diffusion coeffi- 
cient and estimate the diffusion coefficient by its 
value when the mean number of collisions per par- 
ticle is s = 1000, i.e., D — D(s — 1000). 

• We use a smaller system of V = 10 4 particles where 
the long-time tails are cut off due to finite size ef- 
fects. Hence the long-time limit of the integral over 
the VACF is constant and allows to extract the dif- 
fusion constant [30]. 

The diffusion coefficients extracted by these two different 
methods are depicted in Fig. 10. First, we note that the 
results of both methods agree quantitatively, so that it is 
indeed possible to estimate the diffusion coefficient from 
our data. Second, the diffusion coefficient decreases by 
several orders of magnitude as the glass transition is ap- 
proached. Third, for more inelastic systems the decrease 
of D is shifted to higher densities. 

D. Incoherent structure function 

Complete information about the motion of a tagged 
particle is contained in the incoherent scattering func- 
tion, defined as 

1 N 

^c(fc,t) = ^E( eIk - (ri(t) " ri(0)) >- ( 16 ) 

i=l 

In a molecular as well as in a colloidal glass, this func- 
tion displays a two-step decay: The so called /3-relaxation 
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FIG. 10: (Color online) Diffusion coefficients for small parti- 
cles extracted from the different data sets; for detailed expla- 
nation see text. 
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FIG. 12: (Color online) Si uc (k = 3, t) (full lines) for small par- 
ticles in comparison to the Gaussian approximation (crosses); 
parameters as in Fig. 11 
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FIG. 11: (Color online) Si nc (k = 0.5, t) for small particles for 
e = 0.90 and different packing fractions r\ = 0.1; 0.2; 0.3; 
0.4; 0.5; 0.6; 0.65; 0.7; 0.72; 0.74; 0.75; 0.76; 0.77; 0.78; 0.79; 
0.80; 0.81. Inset: Relaxation time r at which the incoherent 
scattering function has decayed to half its initial value as a 
function of packing fraction. 



The two-step relaxation is related to the plateau in 
the MSD and hence more easily observed for larger wave 
numbers. We plot in Fig. 12 the incoherent scatter- 
ing functions S- mc (k — 3, t) for various densities; a two- 
step relaxation is clearly observed for packing fractions 
f] > 0.78. We also compare Si nc (k,t) to the Gaussian 
approximation 

S^ USS (k,t) = e -fc 2 <(Ar(t)) 2 >/4^ (17) 

As expected the Gaussian approximation works well for 
small densities at all times, whereas for the largest den- 
sities only the short time behaviour is approximately de- 
scribed by a Gaussian. The data for the incoherent scat- 
tering function look rather similar to those of the three- 
dimensional system which have been discussed in detail 
in [31]. 



towards a plateau and the subsequent a-relaxation for 
asymptotically long times. The timescale of the a- 
relaxation diverges as the glass transition is approached. 

In Fig. 11, we plot the incoherent scattering func- 
tion for e = 0.9 and increasing packing fraction up to 
rj = 0.81. One clearly observes a pronounced slow- 
ing down which is quantified in the inset of Fig. 11. 
There we plot the relaxation time t(t]), the time at which 
Si nc (k = 0.5, t) has decayed to half its initial value. This 
relaxation time is seen to diverge as the glass transition 
is approached. 



E. Dynamic heterogeneity 

We expect the dynamics to become increasingly het- 
erogeneous as the glass transition is approached. To 
demonstrate this point, we plot in Fig. 13 the MSD of 
the 10% slowest and 10% fastest particles as compared 
to an average over all particles in a system of packing 
fraction 77 = 0.8 and e = 0.9. Whereas the slow particles 
are completely immobile, the fast ones are seen to move 
diffusively covering distances which correspond to several 
radii. We plan to analyze dynamic heterogeneities in a 
future project. 
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FIG. 13: (Color online) MSD for fast and slow particles, either 
restricted to small particles only (dashed) or for all particles 
(full line); the two upper curves refer to the fast particles, the 
two lower curves to the slow ones, and the two middle curves 
to all particles. The MSD at t = 10 3 has been used to identify 
slow and fast particles. 



F. Pair correlation and structure factor 

In this section we want to analyze the structure of a 
dense granular fluid and in particular exclude long range 
crystalline order. Translational order is detected by the 
pair correlation function, g(r), which gives the probabil- 
ity to find a particle in a given distance r, from another 
particle at the origin: 




i) 



(18) 



where hq is number density. For a homogeneous system, 
g(r) approaches a constant as r — > oo which is unity with 
the chosen normalization. Practically, g(r) is calculated 
as follows: 



9(r) 




2nrArnoN 



(19) 

The bin size Ar = i/800 has been chosen such that 
the number of bins and the average number of particles 
per bin remain constant for all packing fractions. This 
corresponds to Ar = 0.83 for the most dilute (ry = 0.1) 
and Ar = 0.30 for the densest (77 = 0.81) system. 

For binary mixtures with big and small particles, there 
are three different pair correlation functions: between 
big particles only, between big and small particles, and 
between small particles only. All three radial distribu- 
tion functions vanish for r/dij < 1, because the particles 
are hard, and have a global peak at the contact point, 



r/di 



1. In Figs. 14 and 15 we show pair corre- 



FIG. 14: (Color online) Pair correlation function for small 
particles for e = 0.90 and different packing fractions, where 
das indicates the sum of two small particles' radii. 
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lations for small and big particles with size ratio 1.25 



FIG. 15: (Color online) Pair correlation function for big par- 
ticles for e = 0.90 and different packing fractions, where dbb 
indicates the sum of two big particles' radii. 



as a function of packing fraction. The data have been 
smoothed, implying a continuous growth of the correla- 
tions for small distances, whereas the raw data strictly 
vanish for r/dij < 1- Increasing the packing fraction, the 
peaks in g(r) get higher and sharper. Simultaneously 
oscillations develop for larger r. At packing fractions 
around rj ~ 0.55, the second peak splits into two sepa- 
rate peaks, a phenomenon that is more pronounced for 
the big particles. This splitting of the second peak near 
r/d^ — 2 is a universal phenomena in glasses, but not a 
good indicator of the transition itself. Our results are in 
good agreement with the experimental ones of ref. [6]. 
We have chosen to study a bidisperse system in order 
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FIG. 16: (Color online) Structure factor for s = 0.90 and 
r\ = 0.81, for small particles. The contour plot of the structure 
factor shows a circular pattern which is characteristic of an 
amorphous system. 



to avoid crystallization. Hence we need to verify that our 
system does indeed remain structurally disordered. We 
do so by mimicking a scattering experiment: We compute 
the static structure factor 

S(k) = (pk(t)p-k(*)) , (20) 

where k is the wave vector and p(k) is the collective den- 
sity variable, defined as 

1 N 

<°xW = -7^E eik ' r3(t) - (21) 

The computed structure factor is shown in Fig. 16 for 
the case of e = 0.90 and r\ = 0.80. The pattern is clearly 
isotropic, revealing the prefered distances which are visi- 
ble as peaks in g(r). We conclude that crytallization has 
been avoided and our system remains structurally dis- 
ordered even when the particles - or a large fraction of 
them - are arrested. 

An alternative measure {Q&) for orientational order 
was introduced in ref. [32] and applied to granular me- 
dia in ref. [33]. We have computed Qq for a range of 
packing fractions and find no indication of long range 
orientational order in agreement with ref. [33] . 

IV. CONCLUSION AND OUTLOOK 

Simulations of a two-dimensional granular fluid in a 
stationary state reveal strong signatures of a glass tran- 



sition, as the packing fraction is increased to approxi- 
mately rj = 0.8. The relaxation time of the incoherent 
scattering function diverges; the MSD develops a plateau; 
the diffusion constant goes to zero; the dynamics becomes 
increasingly heterogeneous. The results of event driven 
simulations are in qualitative agreement with recent ex- 
periments using vibro-fluidized [5] and air fluidized [6] 
granular systems. 

We not only confirm the experimental results, but in 
addition investigate, how glassy dynamics depends on the 
degree of energy dissipation. In particular, we have tested 
the predictions of a recent mode-coupling theory [4] for 
the dependence of the glass transition on inelasticity. The 
theory predicts that the transition is pushed to higher 
densities the more inelastic the system is. To check this 
point, we have simulated a range of coefficients of resti- 
tution, 0.5 < e < 0.9, and compared the mean square 
displacements as a function of e. The time r(e, 77), where 
(Ar 2 (r)) = 1 is shifted to higher densities - in agreement 
with the theory. Similar behaviour is found for the lines 
of equal diffusivities. In addition to the shift in density, 
mode-coupling theory predicts that the plateau in the 
MSD is less pronounced for the more inelastic system. 
This prediction is also born out by our simulations. 

Our simulations clearly show that the diffusion coeffi- 
cient does not exist in a two-dimensional granular fluid; 
an appropriately defined time dependent diffusion coeffi- 
cient grows indefinitely for long times. However, we are 
not able to extract the functional form unambiguously. 
The singular contribution is suppressed near the glass 
transition due to the divergence of the static shear vis- 
cosity, so that our estimates of the diffusion coefficient 
become increasingly more reliable as the glass transition 
is approached. 

Dynamic heterogeneities have been investigated exten- 
sively in molecular glasses [34]. Various tools, such as 
multipoint correlations, have been used to detect increas- 
ing spatial correlations in a disordered structure. For the 
future, we plan to use these methods in order to analyze 
dynamic heterogeneities in a granular medium. 
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